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Abstract. Travelling solitary waves in the one-dimensional discrete nonlinear 
Schrodinger equation (DNLSE) with saturable onsite nonlinearity are studied. A 
variational approximation (VA) for the solitary waves is derived in an analytical form. 
The stability is also studied by means of the VA, demonstrating that the solitons are 
stable, which is consistent with previously published results. Then, the VA is applied to 
predict parameters of travelling solitons with non-oscillatory tails (embedded solitons, 
ESs). Two-soliton bound states are considered too. The separation distance between 
the solitons forming the bound state is derived by means of the VA. A numerical scheme 
based on the discretization of the equation in the moving coordinate frame is derived 
and implemented using the Newton-Raphson method. In general, a good agreement 
between the analytical and numerical results is obtained. In particular, we demonstrate 
the relevance of the analytical prediction of characteristics of the embedded solitons. 
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1. Introduction 

One of the major issues in studies of spatially discrete systems is whether such systems 
can support solitary waves that travel without losing energy to radiation, which results in 
deceleration and eventually pinning of the solitons. The celebrated Peierls-Nabarro (PN) 
barrier pQ is the reason that discrete systems do not generically support exponentially 
localized travelling solitary waves. The barrier corresponds to the energy difference 
between solutions for the on-site- and the inter-site-centred lattice solitons, with the 
latter usually having a higher energy. 

In the class of discrete systems, a ubiquitous model with profoundly important 
applications in physics and applied mathematics is the discrete nonlinear Schrodinger 
equation (DNLSE) [2j. One- and two-dimensional (ID and 2D) equations of this type are 
fundamental models of discrete nonlinear optics, representing planar and bulk arrays of 
nonlinear waveguides coupled by the tunnelling of light between adjacent guiding cores 
[3]. Another well-known application of the DNLSE in any dimension is the description 
of Bose-Einstein condensates trapped in deep optical-lattice potentials, which split the 
condensate into an array of droplets [4] . 

The first attempt at finding travelling lattice solitons in DNLSE was undertaken 
in [5], followed by a more systematic study in [61 [TJ. The latter works indicate that 
travelling lattice solitons in the DNLSE are accompanied by nonzero-radiation tails, 
which was confirmed by more recent studies [H [9j [10] . 

While the most typical onsite nonlinearity in DNLSE models is cubic, waveguiding 
arrays made of photorefractive materials feature the saturable nonlinearity, which 
strongly facilitates the creation of diverse discrete solitons [11] , including solitary vortices 
|12j . necklace-shaped sets p3], circular solitons [2], and symmetry-breaking modes [T5] . 
Various higher-order soliton patterns, in the form of stable multi-charged vortices and 
"supervortices" [16] have also been predicted, in the framework of the same setting. 
Under appropriate conditions, the saturable nonlinearity may be approximated by a 
cubic-quintic truncation; ID and 2D solitons in the DNLSE with the cubic-quintic 
onsite nonlinearity have also been studied in some detail [17]. It was thus found that 
the saturable nonlinearity readily supports travelling solitons in discrete media [18j [19] . 
A reason for this property is the fact that the PN barrier can change its sign in the case 
of the saturable nonlinearity [19], hence the barrier may vanish at isolated points. This 
property may be essential in finding lattice solitary waves that can travel permanently 
without emitting radiation (lattice phonons). 

The saturable DNLSE modelling the propagation of optical waves in a 
photorefractive medium is 

■ dUn a u\ a u\ , au n{t) m 

where u n is a complex-valued wave function at site n, e the strength of the coupling 
between adjacent sites, A 2 u n (t) = u n+ \{t) — 2u n {t)+u n -i{t) the ID discrete Laplacian, A 
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a background frequency, and a the nonlinearity coefficient, which we scale to be a = +1, 
implying that the onsite nonlinearity is self-focusing. To study travelling-wave solutions 
of Eq. (CQ) , an ansatz of the form 

u n (t) = rl>(z,r)e ikn , (2) 

with z = n — ct and r = t, is substituted to yield the time-dependent advance-delay- 
differential equation 

iip T (z, r) = icif} z (z, r) + (2s - A)i/;(z, r) 

- e [Hz + 1, r)e* + - 1, r)e^] + ^^^^ (3) 

where c and are, respectively, the velocity and the wavenumber. 

Travelling-wave solutions of Eq. (CQ) can be sought using the time-independent 
version of Eq. (J3J) , 

ictf + (2e - Aty(z) - £ + l)e ifc + ^(z - l)e~ ik ] + 1+ ^ )[2 = 0> ( 4 ) 

where with ip' = Equation (j3J) takes a simpler form in the case of = 0, following 
Ref. [201 HQ: 

zq// + (2e - Aty(*) - ^ foK* + 1) + 1>{z ~ 1)] + : ^ M2 = °" ( 5 ) 

i + iVvJr 

The existence of travelling solitons for the DNLSE of this type was investigated 
numerically by Melvin et al. j2(Jl EI], using a pseudo-spectral method to numerically 
solve Eq. (JSJ), which yielded weakly delocalized solitary waves. The derealization means 
that the travelling solitary waves in the saturable DNLSE are, in general, accompanied 
by a nonzero oscillating tail, as frequency A will always resonate with the system's 
linear spectral (phonon) band. Because of this, a genuinely travelling solitary wave 
should be an embedded soliton (ES), which can exist inside the continuous spectrum, as 
an exceptional solution (note that the concept of ESs, which was originally developed for 
quiescent solitons in continuous media [22], was later extended to moving pulses [23] and 
to solitons in dynamical lattices [21]). Genuinely localized pulse-like solutions were then 
generated by finding zeros of the amplitude of the soliton's tail (the "tail condition" , 
which is similar to that used to single out ESs in the continuous family of delocalized 
intra-band quasi-solitons [22]). The stability of the numerically obtained solutions can 
then be analyzed in a numerical form by calculating the Floquet multipliers of the 
solutions, using methods similar to that developed in [8]. 

The presence of genuinely travelling lattice solitons in the saturable DNLSE in the 
strong-coupling case has been shown analytically by Oxtoby and Barashenkov, using 
exponential asymptotic methods [25J (see also [10J ) . The use of this sophisticated 
technique is necessary, as the radiation emitted by moving solitary waves is exponentially 
small in the wave's amplitude. This is a reason why broad, small- amplitude pulses are 
highly mobile, seeming like freely travelling solitons. 
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In this paper we apply, for the first time, a variational approximation (VA) to the 
study of travelling solitary waves and their stability, as well as for predicting the location 
of the genuinely localized travelling solitary waves. In the context of DNLSE with cubic 
nonlinearity, the application of VA was proposed to construct the fundamental onsite 
and intersite soliton solutions by using a trial function containing unknown parameters 
that have to be minimized using the Euler-Lagrange equations [261 [27]. The same VA 
method has been applied recently to the cubic-quintic DNLSE [28] (see also [29] and 
references therein). It was shown that the method is not only excellent in approximating 
the fundamental discrete solitons, but also correctly predicts their stability. 

The rest of the paper is organized as follows. In Section 2, we develop the VA 
for the solitary-wave solutions of the advance-delay equation. In the same section, we 
derive an analytical function whose zeros correspond to the location of ESs. The use of 
the VA in analyzing the stability of the travelling solitary waves is then discussed. In 
addition to single-hump pulses, in the same section we consider bound states built from 
two solitons, and the use of the VA to predict the distance between them. In Section 3, 
we introduce a numerical scheme for solving Eq. (jSJ) and compare the numerical results 
with the analytical calculations performed in the preceding section. Results of the work 
are summarized in Section 4. 



2. The variational approximation 

2.1. Core soliton solutions 

As suggested by previous work [8], a travelling lattice wave may be considered as 
superpositions of an exponentially localized core and extended background built of finite- 
amplitude plane waves. Here, we first derive the VA for the core. To this end, we recall 
that Eq. ([5]) can be represented in the variational form, 

5L/5i/j*(z) = 0, (6) 

where 5 /Sip* stands for the variational derivative of a functional, the asterisk denotes 
complex conjugation, and the Lagrangian is 

Hoo r • 

(2e - A)|?/f + In (1 + |^ 2 ) + - [P^ - 

- £ - {rm? + 1) + ^{z - 1)] + ipiriz + 1) + p( Z - 1)]>] d Z . (7) 

A suitable trial function, or ansatz, may be chosen as 

</w0) = F(z) exp (ipz) , (8) 

F(z) = Asech{az), (9) 

where A, a, and p are real variational parameters. While this ansatz postulates 
exponential tails of the soliton, the prediction of solitons within the framework of the 
VA does not necessarily mean that the corresponding solitons exist in a rigorous sense, 
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as the actual tail may be non- vanishing at \z\ — > oo. In fact, the prediction of solitons 
by the VA may imply a situation in which the amplitude of the nonvanishing tail is not 
zero, but attains its minimum |30j . 

The next step is to substitute ansatz ([8]) into Lagrangian ([7]), perform the 
integration, and derive the Euler-Lagrange equations, 

OL/dA = dL/da = dL/dp = 0. (10) 

By substituting ansatz ([8]) into the Lagrangian and performing the integration, we obtain 
the following effective Lagrangian, as a function of parameters A, a, and p: 

2A 2 (2e -A-cp) In 2 (VTTA^ + A) + In 2 (VTTA^ - A) 



L 



<-ir 



a 

\2r 



AA 2 e cos(p) 



sinh(a) 

Then, substituting Lagrangian f lTTj) into Eqs. ffTOj) yields the following equations: 



A(2e-A-cp) In (V 1 + A 2 + A) 2Aecos{p) 
a a\/l + A 2 sinh(a) 

2A 2 {2e-A-cp) In 2 (VTTA 2 + A) + In 2 (VTTA 2 - A) 



a 2 a 2 



| 4A 2 e cos(p) cosh(a) = Q 
sinh 2 (a) 



_ c + 2^ = Q _ 
a smn(a) 

This system of algebraic equations for A, a, and p can be solved numerically. 



2.2. Prediction of the VA for embedded solitons 

We now seek a condition for the possible existence of ESs. To do so, we begin by 
considering a delocalized solution of the linearized version of Eq. ([5]), ^bckg(^), which 
represents the non- vanishing background. Then, following Ref . [30] , it can be shown that 
the condition for the possible existence of ESs, i.e., the absence of nonzero backgrounds 
attached to the soliton, is the natural orthogonality relation, 

{WI* (lH „ M ^bc kg W + cc.} dz = 0, (15) 

where cc. stands for the complex conjugate of the preceding expression and the 
variational derivative 5L / 5ip*\^ z ^ =1 p covc ^ is the left-hand side of Eq. (J5J) with ip(z) 
replaced by 4> core (z). In the context of the VA, T/We^) in Eq- (fT5j) should be substituted 
by the (approximate) form corresponding to the soliton. Here, the background function 
is taken as 

^bckg(z) = ^oexp (iXz) , (16) 
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with constant amplitude ipo, while the soliton's core is approximated by ansatz (181) . Next, 
the frequency A of the oscillating background in Eq. (TT61) can be found by substituting 
■^bckg into the linearization of Eq. (jHJ), i.e., 

icj^ + (2e - A + 1) ip(z) - e [^{z + 1) + ij){z - 1)] = 0, (17) 

which yields 

cA+(A-l)-2e(l-cos(A)) = 0. (18) 

By setting ip = a + if3, where a and (3 are non-zero real constants, Eq. (|T5|) yields 

"+00 

(aM + (3Af)dz = 0, (19) 



where we define functions 

M{z) = |(2e - cA - A) F(z) + ; ^ | cos((A -p)*) 

- e [F(z + 1) cos((A - p)z -p) + F{z - 1) cos((A - p)z + p){20) 
M{z) = | (2e -cX-A)F(z) + 1 ^ } sin((p - \)z) 

- e [F(z + 1) sin((p - \)z + p) + F(z - 1) sin((p - A)z - p)](21) 

It is readily checked that A/"(z) is an odd function, while Ai(z) is an even one. 
Therefore, after some manipulations, integral relation ( fT9l) may be cast into the form 

(2e — cA — A) cos((A — 2 cosh(a2) cos((A — 

cosh(a^) cosh(2az) + 1 + 2A 2 

Be cos((A — p)z) cosh(a2;) Ce sin((A — p)z) sinh(a2; 



dz = 0, (22) 

cosh(2az) + cosh(2a) cosh(2az) + cosh(2a) 

with B = 4 cos(p) cosh(a) and C = 4 sin(p) sinh(a). The integrals in the first and the last 
terms are evaluated using formulas 3.981-3 and 3.983-5, while the second and the third 
terms use 3.984-4, from tables of integrals given in book [31]. The calculation yields 



/N cos ( A - p cosh" 1 1 + 2A 2 /2a) , s 

E = (2e - cA - A) - 2e cos A + V lK . , t- — — ^ = 0, 23 

V ; V; cosh ( [cosh" 1 (1 + 2A 2 )] j2) V ' 

provided that a > and A > p. Thus, in the framework of the VA, Eq. (j23p . along with 
the results of the VA for the soliton's core given by Eqs. (I12p - (]14p . and with Eq. (118)) . 
may predict a curve — in particular, in the (A, e) plane — along which the existence of the 
ESs may be expected. 

2.3. The VA-based stability analysis 

Here, we propose to use the VA to study the stability of the core of the travelling 
lattice solitary wave by calculating eigenvalues for modes of small perturbations in the 
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moving coordinate frame, following Ref. [32] • The stability of the background, i.e., the 
modulational (in)stability of the plane lattice waves, was studied in [8]. 
The underlying time- dependent equation with k = simplifies to 

-ictp z (z,r) +#tO,t) = 

(2e - AMz, r) - e ty{z + 1, r) + ^(z - 1, r)] + 1+ ^ ) r)|2 - (24) 
The Lagrangian of Eq. (I2^j) is 



(2e - A)|^| 2 + In (1 + |^| 2 ) + - - 

oo L ^ 

| + 1, r) + ^(2 - 1, r)] + + 1, r) + ^{z - 1, r)]} 



dz. (25) 



Note that Eq. (EMI) is produced by the variation with respect to ^* not of Lagrangian (12"5|) . 
but rather of the corresponding action functional, S — J Ldt. However, for practical 
purposes (the derivation of VA equations), it is enough to calculate Lagrangian (j2"5|) (it 
is not necessary to calculate the action functional explicitly). 

The time- dependent ansatz, generalizing the static one ([8]), is 

Vw(^, t) = A(r)sech [a(r)(z - f (r))] 

x exp fi<f>(r) + ip(r)z + % -C{r) [z - £(r)] 2 ) , (26) 

where all parameters are real functions of time. Additional variational parameters which 
appear here are the coordinate of soliton's centre, £(t), the overall phase, 0(r), and the 
intrinsic chirp, C(r). Substituting ansatz (1261) into Lagrangian ( 1251) and performing the 
integration yields the corresponding effective Lagrangian, 

-2A + eQ(r) + 2 [£(t)p'(t) + ft (t)-cp(t)] , 7r 2 C"(r) 



J efF 



2 



a(r) 12a(r) 



In 2 (Vl + A(r) 2 + A(r)) + In 2 f ^1 + A(r) 2 - A(r) 
+ — ~~~r~\ — -» ( 27 ) 



with primes standing for the derivatives, and 



47r sin ( ^p- ) cos (p(t)) 

Q{r) = A 



sinh (a(r)) sinh 



C(t)tt 
2a(r) 



_ 4a(r) cos(p(r)) (a(r) 2 + tt) cosGp(t))C(t) 2 4 

sinh(a(r)) 6a(r) sinh(a(r)) 1 1 1 J 

The Euler-Lagrange equations for the variational parameters take the form of 
an ODE system, which may be symbolically written in the vectorial form, x = 
[A / (r),a / (r),p , (r),C ,/ (r),0'(r),«e / (r)] T = g(x), and solved numerically. The VA- 
predicted stability analysis is based on the linearization, 



z = x + <5y, 



(29) 
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with infinitesimal 5, and x = [A, a,p, C = 0, = 0, £ = 0] T representing solutions of 
static variational equations f|T2|) -f lT4|) . The substitution of this into the dynamical 
Euler-Lagrange equations and linearization leads to an eigenvalue problem, 

y = Hy, (30) 

with the corresponding stability matrix H. The stability of the stationary solution is 
then determined by the eigenvalues Q of Eq. ( l30l . which must be found in a numerical 
form, the solution being stable if Re(f2) < for all eigenvalues. 



2.4. The effective potential of the soliton- soliton interaction and the formation of 
bound states 

DNLSEs are known to admit bound states of fundamental (single-humped) solitons 
[33| IM| [35] . The present model also supports bound states, in addition to the single- 
hump solitons. In the infinite domain, there are infinitely many different bound states. 
An essential feature of such states is that the distances between the bound solitons are 
not arbitrary. In the weak-interaction limit, this feature may be explained by means of 
the VA, as we indicate below. 

The effective potential for the interaction between two identical solitons separated 
by distance |/|, which is essentially larger than the width of each soliton, and with a phase 
shift between them (0 = 2 — 1; where <f>i^ are the phases at central points of the two 
solitons), can be derived following the general approach elaborated in Refs. [361 [37]. To 
this end, we consider the wave field in the vicinity of one of the two solitons (say, soliton 
No. 1), whose centre is located at z = 0, while the other soliton (No. 2) is located far 
afield, at z — I] thus we set 

ii>(z) = ti> 1 (z)+1> 2 (z). (31) 

Here ip2 is realized as a weak tail of the second soliton (of course, the tail is affected by the 
overlap with soliton No. 1). Then, expression ( 13~TT) is substituted into the Hamiltonian 
(H) corresponding to Lagrangian ([7]). The effective interaction potential is represented 
by the corresponding term in the total Hamiltonian which is linearized in the weak field, 
ip 2 - Thus, the corresponding contribution to the potential is 



U 



12 



SH , . 5H . 



dz, (32) 



5ip(z) 5ip*(z) 

where the variational derivatives are taken at ip2 — (i.e., for ip = ipi, see Eq. ( !3T|) ). 
The integral is formally written over an infinite domain, although it is assumed that it 
will be calculated in a vicinity of the first soliton (see below). 

Because the stationary one-soliton solution is itself found from the equation 

with ip = ipi, it might seem that expression (132!) should identically vanish (the variational 
derivatives corresponding to the single-soliton solution should be zero, according to 
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Eq. f !33p ). However, the integral will in fact vanish only after performing the integration 
by parts of the terms which contain the first derivative of ip2- 



ic 
~2 



1 dz I dz 



dz. (34) 



Thus, the sole nonzero contribution to U12 originates from the surface term produced 
by the integration by parts in expression (I34I) : U12 = {ic/2) A {ip{ip2 ~ ^1^2} 1 where 
the notation A {...} denotes the difference between the values of this expression at two 
arbitrary points, Z_ and Z + , located sufficiently far from the centre of the first soliton 
(on its left and right sides, respectively), but so that the second soliton remains much 
further still. In fact, one can take Z_ = —00, while Z + is an arbitrary intermediate 
point between the two solitons, located far from both, but closer to the first soliton. 
Finally, the total interaction potential also includes a symmetric contribution from the 
vicinity of the second soliton: 

U int = U 12 + U21 = ^icA ~ M2} + fy&i ~ V^i*}) • (35) 

The main trick (which was employed in Refs. [3HJ EZ]) is to use the asymptotic 
expressions for the wave fields of both solitons taken far from their respective centres 
(i.e., their tails). To this end, we first consider the linearized equation ( |T7|) . with soliton 
tails sought in the form of 

$( z ) = ^°) e Aap2 , (36) 

where constant ip^ is determined by the full nonlinear solution. In fact, expression fl36l ) 
can be viewed as a limiting form of ansatz ([S]) as \z\ — > 00, such that ip^ = 2A and 
A ap = a + ip. Note the similarity between expressions ( l36|) and ( fl~6l) . the only difference 
being that A was real, whereas A ap may be complex. Replacing A by — i\ ap in ( |T8i) . we 
find that A ap then satisfies equation 

icA ap + (2e - A + 1) - 2e cosh(A ap ) = 0. (37) 

It is evident that if A ap is a complex root of Eq. ( |37|) . then — A* p is also a root. Thus, 
the pair of the complex roots may be defined through their real and imaginary parts 
as ±A r + i\i, where A r is chosen to be positive, by definition. In this regard, the 
transcendental complex equation fl37|) may be cast into an explicit form if A r and Aj are 
treated as free parameters, while c and A are considered as unknowns. This approach 
yields 

2e 

sinh (A r ) sin (A.;) , (38) 



A, 

A = 1 + 2e 



cosh (A r ) cos (Aj) + y~ srn h (K ) srn (K) 



(39) 



Aj 

Then Eq. (13^1) yields the tails of the two solitons in the form of 

^i(s) pa ^ m exp (-Ar|«| + iXtz) , (40) 
ip 2 (z) pa ^ (0) exp (-X r \z - l\ + iX { (z - I)) (41) 
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(recall that the centres of solitons No. 1 and 2 are assumed to be located at z = and 
z — I, respectively). The substitution of expression (j4"U|) - (j4Tj) into Eq. fl35|) yields an 
explicit result, which does not depend on arbitrary intermediate point Z + appearing 
in the expression for {7 12 ( nor does it depend on the counterpart of Z + arising in L^i), 
because contributions from Z + cancel out in the final expression (cf. Refs. [361 [37]). We 
thus obtain the following expression for the interaction potential fl35|) : 

U- mt (l) = 2c^ (o) | 2 exp(-A r /)sin(A i /-0), (42) 

where the above-mentioned phase shift is taken into account. In this expression, 
everything is known, in principle (recall that i/j^ is to be found from the full nonlinear 
solution for one soliton), if phase shift is considered as a given frozen constant. 

Strictly speaking, the exponentially decaying potential ( l4"2l) is valid for solitons 
whose waveforms are fully localized. If the actual shape of the solitons features small 
nonvanishing oscillating tails, the asymptotic form of the potential at large values of I 
will change accordingly. 

It is straightforward to see that potential ( 142]) gives rise to a set of local minima 
and maxima (as a function of I), which may correspond to a series of two-soliton bound 
states, as well as to more complex multi-soliton bound states. The extrema of the 
potential are located at points 

l n = ^- arctan f -^J + - + ^ Kn , n ■ sign {\) = 0, 1, 2, 3, . . . . (43) 

The extrema are potential minima for even or odd integers (i.e., for n ■ sign (A $) = 
0,2,4,... or n-sign(Aj) = 1,3,5,..., severally), at c\/\ r < and cAj/A r > 0, 
respectively. Note that the separation between the potential minima, Al = 7r/|Aj|, 
does not depend on the frozen phase shift, <fi. 



3. Numerical scheme and comparisons with analytical results 

To solve Eq. (jSJ) numerically, we use a scheme based on the discretization of the equation, 
resulting in a system of difference equations. We employ central finite differences, 
so that the corresponding Jacobian matrix is sparse. The difference equations are 
then solved using the Newton-Raphson method. This is different from the previously 
used pseudo-spectral collocation method [21], in which the dependent variable if> was 
represented as a Fourier series, whose coefficients were then determined by solving a 
system of algebraic equations obtained by requiring the series approximation to satisfy 
the governing equation at collocation points. 

In the framework of the finite-difference method, with grid size Az,we approximate 
ifj{z) on a finite interval, [—L,+L], as follows: if>(z) — > if>(nAz) = if> n , ip(z ± 1) — > 
if>((n±l/ Az)Az) = if} n ±i/Az- F° r ^'{z) = we use either the central two-point stencil, 
(■071+1 — ^> n _i)/(2Az), or the spectral collocation method [381 [39]. in the following, we 
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describe details of the numerical scheme for the two-point stencil, although the spectral 
collocation method has been implemented too. 

Substituting the above discretizations into Eq. flH]) yields 

-(lpn+1 - 4>n-l) + (2e - A)^n - e(VV»+l/A* + ^n-l/A*) + ; , H i 2 = 0. (44) 



2Az vr ' t+1 r "" v v /rK vr " tl/uz 1 1 l + |<0 n | 

The number of grid points is iV = 2L/Az + 1. We use periodic boundary conditions, so 
that tpN-i+j = i>j and = 4>N-j for j — 1,2, ... , M, where M = 1/ Az. 

Next, we solve the resulting system of nonlinear equations ( |4*4"|) numerically for 
a fixed set of parameters (c, e, A), using the Newton-Raphson method with an error 
tolerance of order 10~ 15 . To do so, we define the left-hand side of Eq. (144)) as f n and 
then seek a solution with f n = for n = 1,2, N — 1. Because ip is complex, we seek 
solutions in the form of ip — Re(ip) + ilm(tp). Accordingly, we define a (real) functional 
vector, F = [Re(/i), . . . , Re(/jv-i), Im(/i), . . . , Im(/jv-i)] T , and a (real) solution vector, 
^ = [Re(^i), . . . , Re(^jv-i), Im(-^i), . . . , Im(^Ar_!)] T . Note that Eq. (jSD has rotational 
and translational invariance. Therefore, to ensure the uniqueness of solutions, we impose 
two constraints, 

/jv±i+i = Re f^-i) - Re (^^±l+i) = 0, ^ 
fE±l +N = Im (ipN±lj = 0. 



These constraints significantly improve the convergence of the Newton-Raphson scheme. 

Strictly speaking, we do not have a rigorous proof of the convergence of the 
numerical scheme outlined above. Therefore, to check the validity of our findings, we 
benchmarked the results against those reported in [21] for the same parameter values. 
The outputs (not shown here) were indeed identical for Az small enough and L large 
enough. Below, we display results for L = 50 and Az = 0.2, which is sufficient for clear 
presentation. We have also used smaller values of Az and larger L, which confirmed the 
robustness of the results. 

3.1. The soliton's core 

For given parameters e, c, and A, we solved variational equations (|T2|) - (fT4|) to produce 
suitable real solutions A, a and p, which then give a quasi-analytical approximation for 
the soliton's core described by function VWe(<z) (see Eq. (JHJ)). As generic examples, 
in Fig. [1] we present, at e — 1 and c = 0.7, the comparison of two soliton profiles 
obtained from the numerical results and VA for two different values of A. We have found 



A w 1.228, a w 0.554, p w 0.377 for A = 0.5 (Fig. 1(a)), and A w 0.685, a w 0.414, 



p ~ 0.368 for A = 0.7 (Fig. 1(b)). Particularly good agreement is observed in both 
cases. 

To further confirm the agreement for different values of A, e and c, in the 
following we compare parameters A, a and p produced by the VA with their numerical 
counterparts. To calculate the latter, we make use of the following relations, which are 
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-10 -5 5 10 -10 -5 5 10 



z z 
(a) e = 1, c = 0.7, and A = 0.5 (b) e = 1, c = 0.7, and A = 0.7 

Figure 1. The comparison of two soliton profiles for two different values of A, as 
indicated in the caption to each panel. The solid lines correspond to numerical results, 
i.e., solutions of Eq. (|44j) imposed with constraints (|45j) . with Az = 0.2 and L = 50 
(the z-axis is truncated to focus the picture on the soliton's core), both the real and 
imaginary parts being shown. The dotted lines are predictions of the variational 
approximation obtained through solving Eqs. (112l) - (|14[) . 



generated by ansatz (jSJ), (jSJ) at z = 0: 

ifj(0) = A, ifj'(0) = ipA, tff"(0) = -A(a 2 +p 2 ), 



(46) 



and take the left-hand sides of these relations from numerical data. Thus, using the 
central finite differences, we obtain the numerical counterparts of A, p, and a: 

Amm = i^E+l , (47) 



Im ( iJ) n+i 



Im ( ipN+i 
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(49) 



The comparison of parameters (A, a, p) obtained numerically (solid line) and from 
the VA (dashed line) is shown in Fig. 2(a) for varying 1/e and fixed (c, A) = (0.7,0.5). 
We observe that the solid and dashed curves are generally close for all the three 
parameters. Nevertheless, we also obtain isolated values of 1/e, which behave as 
singular points. Near the singularities the numerical results deviate very rapidly from 
the predictions of the VA. In fact, at these singular points the numerically obtained 
solutions are strongly delocalized due to the resonance of the oscillating tails with the 
finite size of the computational domains, hence the positions of such singularities depend 
on L, and they may be considered as artifacts of approximating the infinite region by 
the finite domain. 

Similarly to Fig. [TJ where a better approximation is obtained for larger A, we also 



observe in Fig. 2(b) that the variational and numerical curves for A, a, p are closer for 
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A = 0.7 than those in 2(a) In the latter case, the singularities are present too, even 
though they are less pronounced here. 

Thus we can conclude that the VA provides a reliable approximation for description 
of the soliton's core. 



3.2. Embedded solitons 

As shown in Ref. [21], in the general case the numerically obtained solitary waves are 
weakly delocalized, i.e., oscillatory tails with a nonvanishing amplitude are attached to 
them. Nevertheless, as suggested in Refs. [201 EI], it is possible to find solutions with 
vanishing tails (i.e., genuine solitons) by considering quantity 

A, = Im(^(A0), (50) 

which is a signed measure, whose zeros correspond to solutions with vanishing tails. 
In the present case, the periodic boundary conditions lead to Im(-0(A^)) = 0, since the 
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Figure 3. (a) The signed measure A r (cf. Eq. ([51]) ) as a function of 1/e for c = 0.7, 
A = 0.5, Az = 0.2, and L = 50, where A r is zero at e ps 0.737,0.988,1.306 (or at 
1/e ps 1.357, 1.012,0.766), as shown by empty circles, (b) E versus 1/e (cf. Eq. (25|», 
for the same parameter values (A, c) as in panel (a) and for the corresponding solutions 
for (A, a,p) and root(s) A obtained from the Eqs. (|12p - (fTl|) and Eq. (|18p . respectively. 
It is clearly seen that E ^ in the observed domain of 1/e. However, we can conjecture 
that ESs are located near the maximum of E, i.e., at e ps 0.853 (or at 1/e ps 1.172), as 
indicated by the star. 



imaginary part of ^ is odd. Therefore, we modify the signed measure ( 150]) and redefine 
it as 

A r = Re(<ip(N)). (51) 



3( 


a) 


2 





0.737,0.988,1.306 



(i.e., at 1/e ~ 1.357,1.012,0.766). For the parameter values in Fig. 2(b), we show 
the corresponding A r in Fig. 4(a) , from which we conclude that A r vanishes at 
e « 0.828,1.158 (i.e., at 1/e « 1.208,0.864). It is worthy of note that the plot of 
A r also has singularities which occur at exactly the same points as those in Fig. [2] 
obtained from the soliton's core. 

To predict the location of the ESs, we substitute the solution of Eqs. (!T2l-(fl4"l) and 
the root(s) A of Eq. ( JT8|) into Eqs. ([23]) to find E as a function of e, c, and A. Therefore, 
the existence of the ESs can be predicted by seeking for values of the parameters at 
which E = 0. For the parameter values used in Figs. 3(a) and 4(a), curves for E are 
displayed, respectively, in Figs. 3(b) and 4(b) It is seen that E ^ in all the figures, 
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Figure 4. The same as Fig. [3l but for A = 0.7. Zeros of A r are found at 
e ps 0.828,1.158 (or at 1/e ps 1.2077,0.8636), and the maximum of E occurs at 
e 1.020 (or at 1/e ps 0.980). 



i.e., truly localized solitons cannot be directly predicted by the VA. 

However, we can propose a conjecture, based on a "phenomenological" 
consideration of the figures, that there are two zeros of A r on the left and right of 
a maximum of E. For example, in Fig. [3]we have the maximum of E at e ~ 0.853 (i.e., 
at 1/e ps 1.172), which is located between two adjacent numerically found zeros of A r . 
The same phenomenon also takes place in Fig. HI where the two zero-crossing points 
lie between maxima of E, i.e., at e ps 1.020 (i.e., at 1/e ps 0.980) in Fig. HI We have 
also computed (not shown here) the signed measure A r and E for other combinations of 
parameter values, where we observed the same pattern. Thus, we conclude that, with 
addition of a constant, function E may be able to predict the location of the ESs. We 
suspect that the missing constant, which amounts to the shift of the plot for E vertically, 
is related to the choice of the ansatz (see, e.g., [10] for different ansatze accounting for 
the oscillating tails). 



3.3. Stability 

To determine the VA-predicted stability of the soliton, we have solved the 

with 



i.e., 



eigenvalue problem ( 1301) . For the soliton shown in Fig. l(a 
x ps [1.228, 0.554, 0.377, 0, 0, 0] T , we obtain the corresponding eigenvalues 
tt ps 0,0,0,0,0.436z,-0.436z. In addition, for the soliton in Fig. fiTbJ 
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i.e., x ~ [0.685, 0.414, 0.368, 0, 0, 0] T , the corresponding eigenvalues are Q ~ 
0, 0, 0, 0, 0.225z, — 0.225z. As the real part of all eigenvalues is zero, we conclude that 
both solitons are stable. These results are in agreement with the numerical findings of 
Ref. [21]. 

3.4. Bound states 

Two examples of bound states consisting of two solitons, found numerically, are 
presented in Fig. [5] each for the same parameter values. 




-200 -150 -100 -50 50 100 150 

Z 

(a) \l\ ps 39.6 




-200 -150 -100 -50 50 100 150 

Z 



(b) \l\ ps 60.87 

Figure 5. Two bound states, for e = 0.7, c = 0.7, A = 0.7, and Az = 0.2, 
found numerically in the interval of [—200,200], for different distances between the 
two humps, \l\, as indicated in each panel. Solid and dashed lines depict Re(V') and 
Im(^), respectively. 

Next, we compare the prediction presented above with the numerical results shown 
in Fig. For the parameter values used in that figure, we have A r ~ 0.441 and 
A; ~ 0.505. Note that, in terms of the above analysis, the soliton centred at z = in 
Fig. [5] is in fact soliton No. 2. Therefore, phase difference <fi in this case is calculated as 
the phase of the soliton on the right minus the phase of its counterpart on the left. From 
the numerical results, the phase difference between the two solitons shown in the top and 
the bottom panels of Fig. [5] is given, respectively, by ~ —2.012 and ~ 2.368, both of 
which correspond to cAj/A r > 0, hence n-sign (Aj) is odd. Using Eq. (143]) . we find that the 
best fit to the numerically computed distance in Fig. [5] i.e., \l\ ~ 39.6 and \l\ ~ 60.87, 
is provided by l 7 w 41.252 and lg ~ 62.367. Considering the fact that Eq. fj43j) is 
based on a simple approximation (Eqs. ( 1361) . ( 140]) and (HIT) ) and the assumption that 
the interacting solitons have decaying oscillating tails, the approximation is in relatively 
good accordance with the results provided by the numerical solution of the full equation. 
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Using the Runge-Kutta method, we have also carried out numerical integration of 
evolution equation (JTJ, with initial conditions taken as in Fig. |5j The resulting evolution 
of the solutions is displayed in Fig. El where we see that both solitons maintain their 
shapes and positions for a relatively long time. 





Figure 6. Numerical evolution of the solutions presented in Fig.O obtained by means 
of a numerical grid with 401 sites. Upper figures in each panel display the motion of 
the solution through the lattice for the first 100 time units. Shown is a spatiotemporal 
contour plot of the absolute value of the solution. Bottom figures in each panel depict 
the initial (open circles) and final (crosses) profiles of the absolute value of the solution 
after 100 time units of the evolution. 



4. Conclusion 



The aim of this work is to develop the semi-analytical approach to seeking travelling 
solitons, based on the application of the VA to the differential-difference form of the 
DNLSE with the saturable nonlinearity, in the moving coordinate frame. The predicted 
shapes of the solitons are in good agreement with the numerical findings. The VA is 
also extended to examine the stability of the travelling solitons, showing that they are 
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stable, which is consistent with previous work [2T]. Further, the VA was developed to 
predict the locations of the exceptional solutions for genuine travelling solitons with 
strictly vanishing tails. Bound states of two solitons were briefly considered too. In the 
latter case, the VA predicts the distance between two solitons forming the bound state. 
In the numerical part of the work, we have made use of the numerical scheme for the 
DNLSE written in the moving coordinate frame, which is an equation of the differential- 
difference type. Using the Newton-Raphson method, we have confirmed the existence 
of exceptional solutions for travelling discrete solitons (which are "embedded solitons" , 
in this sense), earlier predicted by means of a different numerical algorithm. We have 
compared the analytical results based on the VA and the numerical findings, concluding 
that they are in good agreement. 
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